function showAnimateOfSolution(obj, varargin)
%SHOWANIMATEOFSOLUTION Show the motion of M solution
    
    warning off;
    tList = obj.getTimeSeriesSignal('solution');
    warning on;
    tList = tList.tList;
    %% input check
    p = inputParser;
    % The range of tList, format: [t0, t1]
    p.addParameter('tRange',[min(tList), max(tList)], @(t)validateattributes(t, {'numeric'},{'real','numel',2, '>=',min(tList), '<=',max(tList)}));
    % Number of points in tList
    p.addParameter('npts', 5e4, @(n)validateattributes(n,{'numeric'},{'integer', '>=', 3, '<=', 5e8})); 
    % The length of history trajectory, in unit of frames.
    p.addParameter('HistoryLen',500, @(n)validateattributes(n,{'numeric'},{'integer', 'nonnegative'}));
    % Step size, for example, the first frame is at t = tList(1), the second frame is at t = tList(1+step)
    p.addParameter('step',10, @(n)validateattributes(n,{'numeric'},{'integer', 'positive'}));    
    % Generate an .mp4 movie file
    p.addParameter('mp4_FileName','', @(name)ischar(name) || isstring(name));  
    
    p.parse(varargin{:});
    
    %% open mp4 file for writing
    if isempty(p.Results.mp4_FileName)
        mp4_flag = false;
    else
        mp4_flag = true;
        videoObj=VideoWriter(p.Results.mp4_FileName,'MPEG-4');
        videoObj.FrameRate=30;
        videoObj.Quality=75;
        open(videoObj);
    end
    
    %% get solution
    sol = obj.getTimeSeriesSignal('solution',linspace(p.Results.tRange(1), p.Results.tRange(2), p.Results.npts));
    Mx = sol.getSignalByName('Mx');
    My = sol.getSignalByName('My');
    Mz = sol.getSignalByName('Mz');
    
   
    %% prepare animation stage
    hfMain = figure('position',[300      50      1150         930]);
    ha3D = axes('position', [0.10    0.35    0.53    0.58], 'box','on');
%     ha3D = subplot(3,3,[1,2,4,5], 'box','on');
    haXY = axes('position', [0.69    0.69    0.25    0.25], 'box','on');
%     haXY = subplot(3,3,3, 'box','on');
    haXZ = axes('position', [0.69    0.38    0.25    0.25], 'box','on');
%     haXZ = subplot(3,3,6, 'box','on');
    ha2D0 = axes('position', [0.10    0.05    0.84    0.20], 'box','on');
%     ha2D0 = subplot(3,3,7:9, 'box','on');
    
    % init. 2D diagram
    [ha2D_old, hf2D] = sol.plot();
    ha2D = copyobj(ha2D_old, hfMain);
    set(ha2D, 'position', ha2D0.Position);
    delete(ha2D0); delete(hf2D); clear ha2D_old hf2D ha2D0;
    hold(ha2D,'on');
    tMarker = plot(Mx.tList(1)*[1,1], get(ha2D,'ylim'), 'k--', 'DisplayName', 'Current Time');
    legend(ha2D,{}, 'position', [0.83 0.25 0.07 0.02]);
    
    
    theta = linspace(0,2*pi,200);
    maxAbsM = 1.1*sqrt(max(Mx.signal.^2 + My.signal.^2 + Mz.signal.^2));
    
    scale = 1;
    % XY
    hold(haXY, 'on');
    currXY = quiver(haXY, 0,0, 0,0, scale, 'r');
    hisXY = patch(haXY, 0,0,0 , 'EdgeColor','m','FaceColor','none','EdgeAlpha','flat','FaceVertexAlphaData',[]);
%     hisXY = plot(haXY, 0,0, 'Color',[1,0.7,0.7]);
    plot(haXY,maxAbsM*cos(theta),maxAbsM*sin(theta),'--k');
    set(haXY, 'xlim', [-maxAbsM,maxAbsM], 'ylim', [-maxAbsM,maxAbsM], 'DataAspectRatioMode', 'manual');
    xlabel(haXY,'Mx');ylabel(haXY,'My');

    % XZ
    hold(haXZ, 'on');
    currXZ = quiver(haXZ, 0,0, 0,0, scale, 'r');
    hisXZ = patch(haXZ, 0,0,0 , 'EdgeColor','m','FaceColor','none','EdgeAlpha','flat','FaceVertexAlphaData',[]);
%     hisXZ = plot(haXZ, 0,0, 'Color',[1,0.7,0.7]);
    plot(haXZ,maxAbsM*cos(theta),maxAbsM*sin(theta),'--k');
    set(haXZ, 'xlim', [-maxAbsM,maxAbsM], 'ylim', [-maxAbsM,maxAbsM], 'DataAspectRatioMode', 'manual');
    xlabel(haXZ,'Mx');ylabel(haXZ,'Mz');

    % 3D
    curr3D = quiver3(ha3D, 0,0,0, 0,0,0, scale, 'r');
    hold(ha3D, 'on');
    his3D = patch(ha3D, 0,0,0,0 , 'EdgeColor','m','FaceColor','none','EdgeAlpha','flat','FaceVertexAlphaData',[]);
%     his3D = plot3(ha3D, 0,0,0, 'Color',[1,0.7,0.7]);
    perp3D = quiver3(ha3D, 0,0,0, 0,0,0, scale, '--ro', 'ShowArrowHead', 'off');
    proj3D = quiver3(ha3D, 0,0,0, 0,0,0, scale, '--r', 'ShowArrowHead', 'off');
    plot3(ha3D,maxAbsM*cos(theta),maxAbsM*sin(theta),theta*0,'--k');
    plot3(ha3D,maxAbsM*cos(theta),theta*0,maxAbsM*sin(theta),'--k');
    plot3(ha3D,theta*0,maxAbsM*cos(theta),maxAbsM*sin(theta),'--k');
    % xyz axis
    quiver3(ha3D, 0,0,0, 0.9*maxAbsM,0,0, scale, 'k');
    quiver3(ha3D, 0,0,0, 0,0.9*maxAbsM,0, scale, 'k');
    quiver3(ha3D, 0,0,0, 0,0,0.9*maxAbsM, scale, 'k');
    % xy plane
    [theta,rho] = meshgrid(theta, linspace(0,maxAbsM,5));
    [x,y] = pol2cart(theta,rho);
    surf(ha3D, x,y,0*x, 'EdgeColor','none', 'FaceColor','k', 'FaceAlpha',0.1);
    
    set(ha3D, 'xlim', [-maxAbsM,maxAbsM], 'ylim', [-maxAbsM,maxAbsM], 'zlim', [-maxAbsM,maxAbsM], 'DataAspectRatioMode', 'manual');
    set(ha3D,'CameraPosition', [2.5961 0.4873 0.4234], 'CameraTarget', [0 0 0]);
    xlabel(ha3D,'Mx');ylabel(ha3D,'My');zlabel(ha3D,'Mz');
    
    
    %% main loop
    step = p.Results.step;
    hisLen = p.Results.HistoryLen * step;
    lenSig = length(sol.tList);
    for iT = 1:step:lenSig
        % update animation
        M = [Mx.signal(iT); My.signal(iT); Mz.signal(iT)];
        
        hisRange = max(1,iT-hisLen+1) : iT;
        % XY
        currXY.UData = M(1);
        currXY.VData = M(2);
        hisXY.XData = Mx.signal(hisRange);
        hisXY.YData = My.signal(hisRange);
        if iT < hisLen + 3*step
            hisXY.FaceVertexAlphaData = [linspace(0,0.8,length(hisRange)-1), 0].';
        end
        
        % XZ
        currXZ.UData = M(1);
        currXZ.VData = M(3);
        hisXZ.XData = Mx.signal(hisRange);
        hisXZ.YData = Mz.signal(hisRange);
        if iT < hisLen + 3*step
            hisXZ.FaceVertexAlphaData = [linspace(0,0.8,length(hisRange)-1), 0].';
        end
        
        % 3D
        curr3D.UData = M(1);
        curr3D.VData = M(2);
        curr3D.WData = M(3);
        his3D.XData = Mx.signal(hisRange);
        his3D.YData = My.signal(hisRange);
        his3D.ZData = Mz.signal(hisRange);
        if iT < hisLen + 3*step
            his3D.FaceVertexAlphaData = [linspace(0,0.8,length(hisRange)-1), 0].';
        end
        perp3D.XData = M(1);
        perp3D.YData = M(2);
        perp3D.WData = M(3);
        proj3D.UData = M(1);
        proj3D.VData = M(2);
        
        title(ha3D, sprintf('t = %.3gs',sol.tList(iT)));
        
        % update 2D diagram
        tMarker.XData = Mx.tList(iT)*[1,1];
        
        % write mp4 movie
        if mp4_flag
            frame1 = getframe(hfMain);
            writeVideo(videoObj,frame1);
        end
        
        pause(0.02);
    end
    
    if mp4_flag
        close(videoObj);
    end
end


